Basic usage

from langesim import Simulator

Define the stifness k(t) and center c(t) of the harmonic potential

def k(t):
    """Stifness of the potential.
    Varies linearly from ki=0.5 to kf=1.0 in tf=5.0 time
    units, then remain constant at kf
    """
    tf = 5.0
    ki = 0.5
    kf = 1.0
    if t <= 0.0:
        return ki
    if t < tf:
        return ki + (kf - ki) * t / tf
    else:
        return kf
def c(t):
    """Center of the harmonic potential.
    Varies linearly from ci = -1.0 to cf = 1.0 in tf time units, then remains constant.
    """
    tf = 5.0
    ci = -1.0
    cf = 1.0
    if t <= 0.0:
        return ci
    if t < tf:
        return ci + (cf - ci) * t / tf
    else:
        return cf

Create the simulator

simulator = Simulator(
    tot_sims = 100_000,
    dt = 0.001,
    tot_steps = 10_000,
    snapshot_step = 100,
    k = k,
    center = c
)

The simulator has not performed any simulations yet:

print(simulator.simulations_performed)
0

Run a simulation with the same parameters given in the construction of the simulator

simulator.run()

One simulation has been performed

print(simulator.simulations_performed)
1

The simulation can be accessed at simulator.simulation[0]

simulator.simulation[0].results
{'times': array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ,
         1.1,  1.2,  1.3,  1.4,  1.5,  1.6,  1.7,  1.8,  1.9,  2. ,  2.1,
         2.2,  2.3,  2.4,  2.5,  2.6,  2.7,  2.8,  2.9,  3. ,  3.1,  3.2,
         3.3,  3.4,  3.5,  3.6,  3.7,  3.8,  3.9,  4. ,  4.1,  4.2,  4.3,
         4.4,  4.5,  4.6,  4.7,  4.8,  4.9,  5. ,  5.1,  5.2,  5.3,  5.4,
         5.5,  5.6,  5.7,  5.8,  5.9,  6. ,  6.1,  6.2,  6.3,  6.4,  6.5,
         6.6,  6.7,  6.8,  6.9,  7. ,  7.1,  7.2,  7.3,  7.4,  7.5,  7.6,
         7.7,  7.8,  7.9,  8. ,  8.1,  8.2,  8.3,  8.4,  8.5,  8.6,  8.7,
         8.8,  8.9,  9. ,  9.1,  9.2,  9.3,  9.4,  9.5,  9.6,  9.7,  9.8,
         9.9, 10. ]),
 'x': array([[-1.39320598, -2.01826781, -1.92234229, ...,  0.57260334,
          0.47775196,  0.35274937],
        [-0.61909625, -0.59499455, -0.63164027, ...,  2.01253144,
          0.98817005,  0.05541641],
        [-2.51809684, -2.48784868, -3.03441066, ...,  0.11257338,
         -0.40405033,  0.1429242 ],
        ...,
        [-0.82107004, -0.51178216, -0.81680132, ...,  0.78387936,
          1.1034507 ,  1.02430127],
        [-0.33465449, -0.46067137, -0.56094613, ...,  0.78900911,
          0.32416553,  0.24886806],
        [-2.08355921, -2.09995708, -2.18740501, ...,  0.52437001,
          0.55692324, -0.24688365]]),
 'power': array([[ 0.        ,  0.27180005,  0.25864001, ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        , -0.06782586, -0.05585131, ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  0.4282953 ,  0.66320787, ...,  0.        ,
          0.        ,  0.        ],
        ...,
        [ 0.        , -0.08141434, -0.02097029, ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        , -0.0894174 , -0.06826445, ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  0.29743996,  0.34384373, ...,  0.        ,
          0.        ,  0.        ]]),
 'work': array([[ 0.        ,  0.01329571,  0.03928595, ...,  0.32540684,
          0.32540684,  0.32540684],
        [ 0.        , -0.00695885, -0.012054  , ..., -0.56530878,
         -0.56530878, -0.56530878],
        [ 0.        ,  0.04366937,  0.09401924, ...,  3.97770679,
          3.97770679,  3.97770679],
        ...,
        [ 0.        , -0.00762188, -0.01437032, ...,  0.66022454,
          0.66022454,  0.66022454],
        [ 0.        , -0.01207446, -0.02011314, ..., -0.54664451,
         -0.54664451, -0.54664451],
        [ 0.        ,  0.02338372,  0.05378687, ...,  2.76370786,
          2.76370786,  2.76370786]]),
 'heat': array([[ 0.00000000e+00,  2.33633892e-01,  1.83280729e-01, ...,
         -2.72725620e-01, -2.27688065e-01, -1.54592886e-01],
        [ 0.00000000e+00,  4.66032094e-03, -2.59856982e-03, ...,
          1.04164682e+00,  5.29106841e-01,  9.75155947e-01],
        [ 0.00000000e+00, -2.45718749e-02,  4.92216685e-01, ...,
         -4.16009829e+00, -3.56818263e+00, -4.18657183e+00],
        ...,
        [ 0.00000000e+00,  5.08472014e-02,  9.13533239e-03, ...,
         -6.44874457e-01, -6.62877499e-01, -6.67933247e-01],
        [ 0.00000000e+00, -3.50177862e-02, -5.70389036e-02, ...,
          4.58231926e-01,  6.64349461e-01,  7.18072942e-01],
        [ 0.00000000e+00,  1.44641835e-02,  7.03300090e-02, ...,
         -2.94412106e+00, -2.95907450e+00, -2.27987359e+00]]),
 'delta_U': array([[ 0.        ,  0.24692961,  0.22256668, ...,  0.05268122,
          0.09771877,  0.17081395],
        [ 0.        , -0.00229853, -0.01465257, ...,  0.47633804,
         -0.03620194,  0.40984716],
        [ 0.        ,  0.0190975 ,  0.58623593, ..., -0.1823915 ,
          0.40952416, -0.20886504],
        ...,
        [ 0.        ,  0.04322532, -0.00523499, ...,  0.01535008,
         -0.00265296, -0.00770871],
        [ 0.        , -0.04709225, -0.07715204, ..., -0.08841258,
          0.11770495,  0.17142843],
        [ 0.        ,  0.03784791,  0.12411688, ..., -0.1804132 ,
         -0.19536663,  0.48383428]]),
 'energy': array([[3.86527358e-02, 2.85582342e-01, 2.61219418e-01, ...,
         9.13339519e-02, 1.36371507e-01, 2.09466686e-01],
        [3.62719164e-02, 3.39733904e-02, 2.16193462e-02, ...,
         5.12609957e-01, 6.99738822e-05, 4.46119079e-01],
        [5.76154506e-01, 5.95252003e-01, 1.16239044e+00, ...,
         3.93763007e-01, 9.85678669e-01, 3.67289462e-01],
        ...,
        [8.00398293e-03, 5.12293034e-02, 2.76899152e-03, ...,
         2.33540660e-02, 5.35102405e-03, 2.95275967e-04],
        [1.10671162e-01, 6.35789147e-02, 3.35191174e-02, ...,
         2.22585781e-02, 2.28376113e-01, 2.82099594e-01],
        [2.93525140e-01, 3.31373047e-01, 4.17642022e-01, ...,
         1.13111943e-01, 9.81585096e-02, 7.77359420e-01]])}

Let’s perform another simulation with a different set of parameters.

simulator.run(
    tot_sims = 10_000,
    dt = 0.01,
    tot_steps = 7_000,
    snapshot_step = 500
)

Now 2 simulations have been performed

print(simulator.simulations_performed)
2

The second simulation can be accessed at simulator.simulation[1]

simulator.simulation[1].results
{'times': array([ 0.,  5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60.,
        65., 70.]),
 'x': array([[-0.87221322, -1.0955517 ,  0.51170871, ..., -0.10011611,
         -0.77325558,  0.47818033],
        [-2.44897706,  1.22207542,  0.84702523, ...,  1.63530661,
          1.34227907,  0.60301154],
        [-2.14848429, -0.06627565, -0.2994279 , ...,  0.27821104,
          0.14949696,  0.35752   ],
        ...,
        [-0.88761793,  2.94156042,  0.05477785, ...,  1.63768262,
          1.1136945 , -0.42792164],
        [-0.85540488, -1.44017708,  0.30362297, ...,  2.37718747,
          1.35456978,  1.81298381],
        [-0.83874356,  1.7709441 ,  0.65478371, ...,  2.85621192,
          1.18424655, -0.24347578]]),
 'power': array([[ 0.        ,  1.0561501 ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        , -0.08707466,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  0.48213174,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        ...,
        [ 0.        , -0.5881639 ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  1.27201877,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        , -0.27915072,  0.        , ...,  0.        ,
          0.        ,  0.        ]]),
 'work': array([[ 0.        ,  0.98683144,  0.98683144, ...,  0.98683144,
          0.98683144,  0.98683144],
        [ 0.        ,  0.45242369,  0.45242369, ...,  0.45242369,
          0.45242369,  0.45242369],
        [ 0.        ,  1.23391027,  1.23391027, ...,  1.23391027,
          1.23391027,  1.23391027],
        ...,
        [ 0.        , -1.56171708, -1.56171708, ..., -1.56171708,
         -1.56171708, -1.56171708],
        [ 0.        ,  1.39418006,  1.39418006, ...,  1.39418006,
          1.39418006,  1.39418006],
        [ 0.        , -0.56298211, -0.56298211, ..., -0.56298211,
         -0.56298211, -0.56298211]]),
 'heat': array([[ 0.        ,  1.20475466, -0.87169961, ..., -0.38578608,
          0.58130387, -0.85476592],
        [ 0.        , -0.95264858, -0.96560668, ..., -0.77550008,
         -0.91872984, -0.8985074 ],
        [ 0.        , -0.99519243, -0.71940787, ..., -1.30317466,
         -1.2019866 , -1.35727404],
        ...,
        [ 0.        ,  3.44338808,  2.0052821 , ...,  1.76187921,
          1.56502287,  2.57803975],
        [ 0.        ,  1.5778251 , -1.15693651, ..., -0.45108433,
         -1.33654713, -1.06893565],
        [ 0.        ,  0.85365861,  0.61606835, ...,  2.27924254,
          0.5734546 ,  1.32959721]]),
 'delta_U': array([[ 0.        ,  2.19158609,  0.11513183, ...,  0.60104536,
          1.5681353 ,  0.13206552],
        [ 0.        , -0.50022488, -0.51318299, ..., -0.32307638,
         -0.46630615, -0.44608371],
        [ 0.        ,  0.23871784,  0.5145024 , ..., -0.06926439,
          0.03192367, -0.12336377],
        ...,
        [ 0.        ,  1.881671  ,  0.44356503, ...,  0.20016213,
          0.00330579,  1.01632267],
        [ 0.        ,  2.97200516,  0.23724355, ...,  0.94309573,
          0.05763293,  0.3252444 ],
        [ 0.        ,  0.29067649,  0.05308623, ...,  1.71626043,
          0.01047248,  0.7666151 ]]),
 'energy': array([[0.00408237, 2.19566846, 0.11921419, ..., 0.60512773, 1.57221767,
         0.13614788],
        [0.52488363, 0.02465875, 0.01170064, ..., 0.20180725, 0.05857748,
         0.07879992],
        [0.32975404, 0.56847188, 0.84425644, ..., 0.26048965, 0.36167771,
         0.20639027],
        ...,
        [0.00315743, 1.88482843, 0.44672246, ..., 0.20331956, 0.00646322,
         1.0194801 ],
        [0.00522694, 2.97723209, 0.24247049, ..., 0.94832266, 0.06285987,
         0.33047134],
        [0.00650091, 0.2971774 , 0.05958714, ..., 1.72276134, 0.01697339,
         0.77311601]])}

Let’s analyse the results of the first simulation

sim = simulator.simulation[0]

Plot trajectories number 1, 2, 5 and 10

sim.plot_sim('x', [1, 2, 5, 10])

One can plot the following quantities:

sim.result_labels
['x', 'power', 'work', 'heat', 'delta_U', 'energy']

Plot the work done of trajectories 0 to 9

sim.plot_sim('work', range(10))

Plot the average position, the variance of the position, the average work and heat.

sim.plot_average('x')
sim.plot_variance('x')
sim.plot_average('work')
sim.plot_average('heat')

Plot an animation of the PDF of the position

sim.animate_pdf('x', x_range = [-5, 5], y_range= [0, 0.5])

Plot the PDF of work

sim.animate_pdf('work', x_range = [-2, 7], y_range = [0, 1.5])

To perform simulations with a different potential, a new simulator has to be created. Let us create one for a non harmonic potential.

def double_well_potential(x,t):
    """Quartic potential with a double well with moving center"""
    tf = 5.0
    v = 0.5
    if t < tf:
        return (x-v*t)**4 - (x-v*t)**2
    else:
        return (x-v*tf)**4 - (x-v*tf)**2

simulator2 = Simulator(
    tot_sims = 100_000,
    dt = 0.001,
    tot_steps = 10_000,
    snapshot_step = 100,
    harmonic_potential= False,
    potential = double_well_potential
)
simulator2.run()
sim2 = simulator2.simulation[0]
sim2.animate_pdf('x', x_range=[-2, 4.5], y_range = [0, 1])

By default, the initial condition is at thermal equilibrium at t=0. But this can be changed. Let’s perform some simulations with a fixed initial condition.

def init_condition():
    """Initial condition particle is at x=4.0 at time t=0"""
    return 4.0

Let’s create a new simulator with that initial condition and the quartic potential

simulator3 = Simulator(
    tot_sims = 100_000,
    dt = 0.001,
    tot_steps = 10_000,
    snapshot_step = 100,
    harmonic_potential= False,
    potential = double_well_potential,
    initial_distribution = init_condition
)
simulator3.run()
simulator3.simulation[0].animate_pdf('x', x_range = [-2, 4.5], y_range = [0, 0.8])